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ABSTRACT 

The stellar ejection rate and the rates of change of the binary semimajor axis and 
eccentricity are derived from scattering experiments for the restricted three-body 
problem. They are used to study the evolution of binaries in simple models for galac- 
tic nuclei, starting soon after the black holes become bound and continuing until the 
evolution is dominated by the emission of gravitational radiation, or until the ejected 
mass is too large for the galaxy to be considered fixed. The eccentricity growth is 
found to be unimportant unless the binary forms with a large eccentricity. The scat- 
tering results are compared with predictions from Chandrasekhar's dynamical-friction 
formula and with previous work on the capture and scattering of comets by planetary 
systems. They suggest that a binary with masses mi > m2 should not be considered 
hard until its orbital velocity exceeds the background velocity dispersion by a factor 
that scales as (1 -|- mi/m2)^/^. 



1 INTRODUCTION 

The existence of massive black hole (BH) binaries follows 
from two widely-accepted assumptions: that galaxies merge 
with other galaxies, and that many galaxies contain mas- 
sive BHs. For if two BHs enter the core of a merged galaxy, 
dynamical friction drags them to the center where they 
form a binary. The subsequent evolution was first outlined 
by Begelman, Blandford, and Rees (1980, hereafter BBR). 
Initially the binary hardens (i.e. its separation shrinks) be- 
cause of the interaction between the BHs and all the stars 
in the galaxy core. But that is ineffective once the BHs be- 
come close because distant stars perturb the binary's cen- 
ter of mass but not its semimajor axis. The binary then 
hardens by giving kinetic energy to stars that pass in its 
immediate vicinity; a hard binary can eject stars out of 
the core at high velocity. If there are enough stars for the 
hardening to continue (and gas accretion onto the BHs can 
help), eventually the BHs merge through the emission of 
gravitational radiation; otherwise the hardening stalls and 
the binary survives for the lifetime of the galaxy. 

Whether a binary merges or survives and how long it 
spends in each stage of the evolution are questions rele- 
vant to a number of problems in extragalactic astronomy. 
Their answers would help us predict the total BH merger 
rate and whether it is high enough for us to detect the 
resulting gravitational waves (e.g. Thorne 1992, Haehnelt 
1994). They would help us assess BH-binary models for the 
bending and apparent precession of radio jets from active 



galactic nuclei, first proposed by BBR. And they would tell 
us what to expect if three or more BHs enter the core of a 
galaxy, which can happen if the BHs are dragged in from 
the galaxy's halo or if the galaxy undergoes multiple merg- 
ers with other galaxies containing BHs. If the first binary 
merges fast it can form a binary with a third BH, and once 
that merges it can form a binary with a fourth, and so on, 
leading to a massive central BH; but if the first binary still 
exists when a third BH enters then one or all three of the 
BHs can be ejected in a sling-shot interaction. Arguments 
like these can set limits on massive BHs as dark-matter 
candidates for galactic halos (see Hut and Rees 1992, Xu 
and Ostriker 1994 for conflicting limits for our Galaxy). 

Another question is what a binary merger does to the 
surrounding galaxy, i.e. what observable signature it leaves. 
Mass ejection during the evolution should reduce a galaxy's 
central density and expand its core (BBR). Ebisuzaki, 
Makino, and Okumura (1991) have proposed this as an ex- 
planation for why large elliptical galaxies have lower central 
densities and weaker density cusps than small ellipticals 
(e.g. Kormendy et al. 1994). 

We are far from having precise answers to any of these 
questions. BBR gave a range of merger times for one typical 
example that spanned three orders of magnitude because of 
the uncertain influence that mass ejection has on the hard- 
ening rate. Fukushige, Ebisuzaki and Makino (1992) have 
argued that dynamical friction causes a binary to become 
highly eccentric and that this greatly reduces the merger 
time because gravitational radiation then becomes impor- 
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tant early in the evolution. Although their arguments are 
not convincing, they have called attention to the eccen- 
tricity growth and our ignorance of its correct description. 
There are uncertainties in how the hardening rate depends 
on the ratio of the two BH masses, in when a binary makes 
the transition from soft to hard, and even in what the words 
soft and hard should mean in this context. And our knowl- 
edge of how a binary merger changes a galaxy is based on 
back-of-the-envelope estimates and simple N-body experi- 
ments with unrealistic galaxy models. 

The first step towards resolving these questions is to 
understand how a massive binary evolves in fixed stellar 
background. Consider a binary with masses mi > m2 and 
semimajor axis a in an isotropic background of stars of 
mass m« -C m2, density p, and one-dimensional velocity 
dispersion cr; let M12 and /U denote the total and reduced 
binary mass: 



M12 



m-i + m2, 



^ = rairn'z/Mr. 



(1) 



The binary evolution and its effect on the galaxy are de- 
termined by three dimensionless quantities: the hardening 
rate 



Gpdt \aJ ' 



(2) 



the mass ejection rate (where A/cj is the stellar mass that 
the binary has ejected from the galaxy core) 



M12 dln(l/a)' 
and the eccentricity growth rate 

de 



K = 



dln(l/a) 



(3) 



(4) 



The quantities H, J, and K can be found from scattering 
experiments that treat the star-binary encounters one at a 
time. Analytic approximations such as the impulse approx- 
imation are helpful during the early stages of the evolution 
(Gould 1991), but not once the binary becomes hard. 

Most published scattering experiments assume the bi- 
naries and stars to have equal or nearly equal masses (see 
Heggie 1988 for a review). There are some exceptions. Roos 
(1981) performed scattering experiments for the restricted 
three-body problem to study the evolution of hard, massive 
BH binaries, and used them to correct a misplaced factor of 
mi /m2 in the BBR hardening rate. He tried to measure K 
but his statistics were too poor to give definite conclusions 
(only 500 orbits per meeisurement) . Hills (1983a) used the 
general three-body problem to study interactions between 
a massive binary and low-mass intruders (m«/m2 = 0.01). 
He gave results for the hardening rate for a wide range of 
mass ratios {m\/m2 = 1 - 300), but like Roos he considered 
only very hard binaries. Mikkola and Valtonen (1992) used 
the restricted three-body problem to measure H and K for 
equal-mass BH binaries with varying degrees of hardness. 
Their measurements are accurate for hard binaries but have 
large error bars for binaries that are not hard. 

If mi ^ m2 then the interaction between a star and a 
BH binary is similar to the interaction between a comet and 
a planet orbiting a star. Although scattering experiments 



are used to study cometaxy dynamics (see Fernandez 1993 
for a review) they are not of much help for our questions 
about BH binaries, partly because they often consider only 
one mass ratio (for the Sun- Jupiter system), but mostly 
because they are used to answer different questions, such 
as the cross section for the capture of interstellar comets, 
or for the conversion of long-period comets to short-period 
comets, the survival probability of comets once they are 
captured, and how all these depend on the comet's inclinar 
tion. There is nevertheless some overlap between the two 
problems. 

The goal of this paper is to present accurate measure- 
ments of H , J, and K over the range of parameters of inter- 
est for the BH-binary problem, including the dependence 
on the mass ratio, eccentricity, and degree of hardness of 
the binary. Other quantities to be studied include the cross 
section for a binary to capture stars into bound orbits, 
for close encounters between stars and the binary mem- 
bers, and the distribution of velocities with which stars 
are expelled from the binary. These add to our under- 
standing of H, J, and K, and are needed by themselves 
for some applications. The results will be presented in a 
model-independent way so that they can be applied to any 
problem with m* <C m2 < mi . 

Once H, J, and K have been measured they can be 
used to study the evolution of binaries in fixed galaxy mod- 
els. That will be done here for some simple models. If the 
BHs are large they will of course eject too much mass from 
the galaxy core for it to be considered fixed. But the results 
will still be valid during the early stages of the evolution. 
And they will be helpful even in the later stages, because 
we can imagine at any instant that the binary is embedded 
in a fixed background whose properties are those of the 
galaxy at that instant. The self-consistent evolution of a 
massive binary in a realistic galaxy model and the changes 
this induces in the model are best studied by large N-body 
experiments. That will be deferred to paper II, along with 
a discussion of what both papers imply for the astronomi- 
cal questions mentioned above (Quinlan and Hernquist, in 
preparation) . 



2 COMPUTATIONAL METHOD 

2.1 Derivation of results from the restricted 
three-body problem 

We treat the star as a massless test particle moving in the 
potential of the two BHs. From the changes in the star's 
energy and angular momentum per unit mass, A_E« and 
AL« , we infer the corresponding changes AE and AL that 
the binaxy would have suffered if the star had been given 
a small but nonzero mass. The three bodies are treated as 
point masses and gravitational radiation is ignored. 

In a real galax;y stars approach the binary with a wide 
distribution of velocities at any given time. But the scat- 
tering experiments are easiest to perform if the stars all 
start from the same velocity v at a large separation from 
the binary (the initial velocity, or the velocity at infinity). 
In that case we write 
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where the "1" indicates that the stars all have a single 
velocity v. The hardening rate for a Maxwellian velocity 
distribution is then 



H{a) = f 
Jo 

where 



2 ^ 

dv Awv f(v, a) — Hi(v) 

V 



1 
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exp(— v^/2cr^). 



(6) 



(7) 



The hardening rate is derived from the average energy 

change for stars that scatter off the binary. We define a 
dimonsionless energy change C by (Hills 1983a) 



C ■ 



Mi2 AE aAE^ 



2m, E 



Gfj. 



(8) 



This must be averaged over all angular variables describing 
the binary's orientation and phase, to give (C), and then 
integrated over all impact parameters. The averaging and 
integrating are done in a Monte Carlo fashion by picking 
orbits from suitable distributions. We sometimes describe 
an orbit by its impact parameter b, the distance at which it 
would pass the binary if it felt no attraction, and sometimes 
by its pericenter distance r^, the distance if it is attracted 
by a point mass M12 . The two are related by gravitational 
focusing: 



1 + 



2GM12 



(9) 



We also define a dimensionless impact parameter x by 

X = b/bo, bl = 2GMi2a/v'^; (10) 

60 is approximately the impact parameter corresponding 
to Tp = o if gravitational focusing is important. With this 
notation we can write 



Hi = 8-Kh{C) = 8n dx x{C) 



f 

Jo 



(11) 



where the second equality defines the operator I^- 

The derivation of the eccentricity growth rate is sim- 
ilar. The change to the binary's eccentricity from a single 
scattering event is, if the change is small and the binary's 
orbital angular momentum points in the z direction, 



Ae 



(1 



(12) 



Aln(l/a) 2e 
We define a dimensionless angular-momentum change B by 

M12 AL^ M12 AL.„ 



m. fj. [GMi2a(l-e2)]i/2' 

and can then write 



(1 



2e 
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where the "1" has the same meaning as before. The deriva- 
tion of K from Ki will be described later. 



2.2 The scattering experiments 

Each scattering experiment requires five uniformly- 
distributed random numbers (four if the binary is circular) : 

one for the square of the impact parameter (in some range 
[0,6max]), and four to fix the binary's orientation and phase: 
the cosine of the inclination ([-1,1]), the longitude of as- 
cending node ([0,27r]), the argument of pericenter ([0,27r]), 
and the mean anomaly at some fixed time ([0,27r]). The 
numbers are chosen with the quasi-random number gener- 
ator sobseq of Press et al. (1992). 

The range [0,6max] for impact parameters is split into 
five intervals corresponding to ranges in scaled pericenter 
distance rp/a of [0,1], [1,2], [2,4], [4,8], and [8,16]. Each out- 
put quantity is measured in a number of steps. On the first 
step the program spends short but equal amounts of cpu 
time picking orbits from the five intervals. On each suc- 
cessive step the program doubles the cpu time and adjusts 
its strategy so that the time it spends on each interval is 
proportional to the uncertainty that interval contributes to 
the quantity being measured. Once the uncertainty is re- 
duced to an acceptable level, or the cpu time exceeds some 
maximum allowed value, the results from all five intervals 
are combined with appropriate weights for a distribution 
uniform in b^. For H, J, and K the last three intervals 
contribute little because the changes in energy and angular 
momentum fall off rapidly with increasing impact parame- 
ter. 

The coordinates are chosen so that the binary's cen- 
ter of mass is at the origin and the star starts at infinity 
with {x,y,z) — (f), 0, 00) and {vx,Vy,Vz) = (0,0, — «). The 
star is moved from r — 00 to r = 50a along a Keplerian 
orbit about a point mass M12 at the origin. The numerical 
integration starts at r = 50a. 

The orbits are integrated in double precision with an 
explicit, embedded Rungo-Kutta method of order (7)8: the 
program dopriS of Hairer, Norsett, and Wanner (1987). 
The program adjusts the integration stepsize to keep the 
fractional error per step in the position and velocity below 
some level e, which was set to 10~®. With this choice the 
change in a star's Jacobi constant for a circular binary is 
at most 10~^GMi2/a and often much smaller. The forces 
from the BHs are not softened. 

Some integrations are time consuming because the star 
gets captured into a weakly-bound orbit and makes many 
revolutions before it is expelled. The integration stepsize is 
a small fraction of the binary's period even if the period 
of the star about the binary is much longer. The following 
approximation expedites those experiments. If a captured 
star moves further than = a(10^''m2/mi)^'''' from the 
binary, the binary is replaced by a point mass M12 and the 
star is moved along a Keplerian ellipse until it returns inside 
ffc, when the forces from the BHs are reintroduced with 
the correct orbital phase. The (m2/mi)^'^'' mass scaling 
makes the quadrupole force from the binary at rk about 
10~^''GMi2/o^, independent of m^/mi. 

Orbits that get captured for long times tend to be 
highly chaotic. The integration for any particular orbit of 
that type is difficult to justify because a small change to 
the integration procedure can make a big change to the out- 
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come. But the average results derived from a large number 
of integrations can be correct even if the individual inte- 
grations are not; that is suggested by shadowing lemmas 
that have been proved for simple chaotic systems. The av- 
erage results presented here do not change noticeably if e 
is raised or lowered by a factor of 100, even though some 
orbits undergo big changes. 

An integration is stopped when the star leaves the 
sphere r — 50a with positive energy. The average results 
are not sensitive to the location of this sphere provided that 
it is at least 10-15 times larger than a. Once the integra- 
tion stops the program records the changes to the star's 
energy and angular momentum, the minimum separations 
between the star and the two BHs, and between the star 
and the binary's center of mass, the integration time, the 
number of integration steps, and the number of times the 
star's radius passed through a minimum (which, if greater 
than one, gives the number of revolutions that a captured 
orbit made) . An orbit that does not get captured typically 
takes a few hundred to a few thousand integration steps. 
Captured orbits can take much longer. If an integration 
lasts for more than 10^ steps it is abandoned. The fraction 
of abandoned integrations is the largest for hard binaries 
with mi/m2 3> 1, but even for those it is less than 0.1%. 

The error in any average quantity has a systematic 
component and a statistical component. Systematic errors 
arise, for example, because of errors in the numerical in- 
tegration, because integrations are abandoned if they take 
too long, because the program imposes a maximum impact 
parameter bmax, and because the orbits start and end at 
r = 50a instead of at r = oo. But none of these is large: the 
total systematic error is usually much smaller than the sta- 
tistical error. The statistical errors are estimated by taking 
the difference between results found with A'^ orbits (the final 
number) and N/2 orbits, or sometimes — if that difference 
looks suspiciously small — one half the difference between 
A'^ orbits and A/4 orbits. That gives a rough estimate of 
the error level. The statistical errors decrease at large A as 
(In Nf/N when quasi-random numbers are used, where d 
is the number of numbers picked for each experiment (5 in 
general, 4 if the binary is circular). That is faster than the 
N~^^^ decrease that occurs with random numbers (Press 
et al. 1992). 

The number of orbits needed to reduce the statistical 
errors to an acceptable level varies widely with the binary 
and the measurement. Measurements are more difficult for 
Ki than for Hi because of cancellation. Cancellation is a 
problem for Hi too at high v values. And regardless of what 
is being measured, binaries with mi ^ 7712 require many 
more orbits than equal-mass binaries because of the rare 
close encounters with the mass m^. When the results are 
presented the number of orbits used will not be given be- 
cause it is different for each measurement; some estimate of 
the statistical error will be given instead. The numbers are 
about 10''-10^ orbits per measurement for Hi (sometimes 
more at high v values) and 10^-10® for Ki. The complete 
set of experiments took about four months of cpu time on 
an IBM 580 RISC workstation. 



mi/m2 


Ho 


«'/Vbin 


1 


17.97 


0.5675 


4 


20.54 


0.4263 


16 


21.87 


0.2228 


64 


22.78 


0.1043 


256 


22.57 


0.0573 



Table 1 . Parameters for fits to Hi (oq. p^ ) for a circular binary. 



3 RESULTS FROM THE SCATTERING 
EXPERIMENTS 

3.1 Hardening 

3.1.1 Hardening rate 

The hardening rate Hi ^ has been measured as a function 
of the binary's eccentricity and hardness for a wide range 
of mass ratios. It is plotted in Figure |l| versus the hardness 
as given by the ratio of the initial stellar velocity v and the 
binary's orbital velocity Vbin (the relative velocity of the 
two BHs if the binary is circular); 



Vbin = ^/GMl2/a. 



(15) 



The error bars show the statistical-error estimates; if not 
visible then they are smaller than the size of the points. 

The velocity dependence of Hi is fit by a function with 
two free parameters (whose values are given in Table ^ : 

Ho 



Hi = 



[1 -f- (v/w) 



411/2 ■ 



(16) 



The function has a constant value Ho at v <^ w, starts 
to decrease as v approaches w, and decreases as l/n^ at 
V ^ w. This fits the data well at high and low veloc- 
ities. It does not fit so well when Hi first starts to de- 
crease; those data points were given little weight in the fit- 
ting procedure. A three-parameter function was also tried. 
Hi = Hq/ [1 + {v/w)''f'''', but the exponent k never dif- 
fered by much from 4. 

For an equal-mass binary the constant hardening rate 
at low velocity. Ho — 18.0, agrees well with Mikkola 
and Valtonen's (1992) nRa — 18.2, and also with the 
resuhs of Hills (1983a, 1992). At high velocity Hi ^ 
Ho{w/vf ~ 5.8(Vbin/^;)^ which can be compared with 
Gould's (1991) analysis using the impulse approximation, 
-ffi, Gould = (87r/3)(Vbin/«)^. The agreement is satisfactory 
considering that the error bars are large at high veloc- 
ity and that the impulse approximation is justified only 
if v/Vbin 3> 1. In fact it is surprising how well the impulse 
approximation works when v/Vhin — 1- 

Panel (a) shows that the hardening rate for an equal- 
mass binary does not vary much with the eccentricity. 
Mikkola and Valtonen (1992) reached the same conclusion, 
which remains true for all mass ratios. For later applica- 
tions the hardening rate for a circular binary will be used 
for all eccentricities because the variation of Ho and w with 
eccentricity is too small to matter. 
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Figure 1. Binary hardening rate versus the initial stellar velocity v. (a) shows three eccentricities for m\/m2 = 1; (b) and (c) show 
five mass ratios for e = 0, on linear and logarithmic scales. The lines are the fits to eq. {[19). 



Roos (1981) and Hills (1983a) showed that the low- 
velocity hardening rate Ho does not vary much with the 
mass ratio. But the velocity w does. The variation is fit by 



0.85 



— yM„ = 0.85W^. 



(17) 



The physical significance of this mass dependence is the 
following: ii v < w the binary can easily capture stars into 
bound orbits; if u > ui it cannot. 

The integral (^ for a Maxwellian distribution was eval- 
uated numerically. The relation between H and Hi is fit 
closely by the formula 



Hi{V3a) 



In 



1 + a 



(18) 



with a — 1.16 and (3 — 2.40. In the limit of high velocity 
this gives 



H : 



f3Ho 
6 



In 



(19) 



The log term looks like a familiar Coulomb logarithm but 
comes from an integral over the velocity distribution, not 
over a range of impact parameters. The limit ( |l9| ) can be 
compared with Gould's (1991) hardening rate, 



16V27r 



bin 



(20) 



For an equal-mass binary the coefficients of the log terms 
differ by about 30%, which is satisfactory considering 
the uncertainties mentioned above. For non-equal masses 
Gould's log term does not have the correct mass depen- 
dence (Gould did not attempt to compute the log term 
accurately) . 

The hardening rate for a massive BH binary is some- 
times derived from Chandrasekhar's dynamical-friction for- 
mula (e.g. Fukushige et al. 1992) . The error in that has been 
known for many years (Chandrasekhar 1944, Hills 1983a): 
the distant encounters included in the friction formula do 
not perturb the binary's semimajor axis — they only per- 
turb its center of mass. It is an accident that the deriva- 
tion gives a result like Gould's for a Maxwellian distribu- 
tion if a suitable choice is made for the log term: if the 
same derivation is used for Hi it gives the nonsense result 
that (for mi — 7712) Hi is zero at w = 0, rises as Hi ~ v 
for V < Vbin/2, and then drops abruptly back to zero at 
V ~ Vbiii/2 (because only stars moving slower than the BHs 
contribute in Chandrasekhar's formula). See Gould (1991) 
for further discussion. 

The velocity dependence of the hardening rate sug- 
gests a new convention for the use of the word hard. A 
hard binary is usually defined in one of three ways. The 
first says that a binary with binding energy is hard if 
-Eb 3> rrifa^ and soft if Ei, <^ mtcr^ (p. 534 of Binney 
and Tremaine 1987). The second says a binary is hard if 
it grows harder through interactions with stars and soft if 
it grows softer (Hut 1983). And the third, which is often 
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stated as a corollary of the first or second rather than as 
an independent definition, says a hard binary is one that 
"hardens at a constant rate," i.e. at a rate H that is inde- 
pendent of the hardness. The equivalence of the first two 
definitions is called Heggie's law. But neither of those defi- 
nitions is useful for a massive BH binary because both are 
satisfied by almost any pair of massive BHs that is close 
enough to be called a binary. 

The third definition gives almost the same result as 
the first two when the masses are equal (see Ifut 1983, or 
Figure 6.3 of Spitzer 1987) and is far more useful when the 
binary is massive. A BH binary will therefore be called hard 
if it hardens at a constant rate, i.e. if u w or equivalently 
if Vbin/o" ^ (1 -I- mi/m2)^^^. It is tempting to call a binary 
soft if it is not hard, but that is confusing for massive bi- 
naries because there is a wide gap for them between a hard 
binary in the sense used here and a soft binary in the famil- 
iar sense that "soft binaries grow softer." In later figures 
the properties of hard binaries are studied with scattering 
experiments using the lowest initial velocity v in Figure hi 
for each mass ratio; those velocities are logj^Q(ii) = —2.025, 
-2.25, -2.6, -2.8, and -3.0 for mi/m2 = 1, 4, 16, 64, and 
256. 

The scattering results and Gould's (1991) analysis re- 
fute Hills's (1990) statement that a binary grows harder if 
Vbin > f and softer if Vbin < o" regardless of the values of 
m,, mi, and m2. Although the mean energy change (C) at 
zero impact parameter does change from positive to neg- 
ative when the stellar velocity is raised from v < Vbin to 
V > Vbin, that sign change disappears when (C) is averaged 
over impact parameter. 

The reason for the hard/not-hard transition at ct = ui 
is best explained after we have examined the cross sec- 
tion for stars to be captured by a binary, to have close 
encounters with the binary members, and the distribution 
of velocities with which stars are expelled from a binary. 



3.1.2 Capture cross section 

We say that a binary captures an incoming star if the 
star's orbital radius passes through more than one mini- 
mum. Almost all captured orbits are eventually expelled in 
the three-body problem (there might be a set of measure 
zero that remain bound forever), but the star can survive 
for many revolutions before that happens. 

Previous work has used scattering experiments and ap- 
proximate methods to derive capture cross sections. Hills 
has used scattering experiments to study the capture of 
orbits by very hard, massive binaries (Hills 1983a, 1983b, 
1992). He unfortunately defines capture — or what he calls 
long-term capture — in a way that depends on his program 
(he says a long-term capture occurs if the integration takes 
more than 150,000 steps). But he gives helpful informa- 
tion on how the capture probability depends on the impact 
parameter, eccentricity, and binary mass ratio. Pineault 
and Duquet (1993) have used the impulse approximation 
to derive approximate capture cross sections for massive, 
circular binaries, for arbitrary mass ratios and degrees of 
hardness (they give many relevant references to the comet 



mi/m2 


Cl 


C2 


C3 


C4 


C5 


1 


17.97 


1.0066 


3.5745 


2.0865 


0.6100 


4 


20.54 


0.7929 


4.5326 


1.2675 


1.2377 


16 


21.87 


0.4122 


3.6588 


1.2324 


0.9754 


64 


22.78 


0.1800 


6.1855 


0.5562 


1.0087 


256 


22.57 


0.0846 


8.1992 


0.3856 


0.9782 



Table 2. Parameters for fits to Ecap (eq. |22[ ) for a circular bi- 
nary. 

literature). They say their cross sections are accurate to 
within a factor of 2-3, although that is not clear because 
they adjust their formulas in an ad hoc way — using Hills's 
(1983a) results to guide them — for hard binaries for which 
the impulse approximation does not work. 

The measurements made here improve upon those of 
Hills by using a reproducible capture definition and by ex- 
ploring the dependence on the binary's degree of hardness. 
Panels (a) and (b) of Figure ^ show the capture cross sec- 
tion for a circular binary in units of the binary's geomet- 
rical cross section Ebin, which includes the correction for 
gravitational focusing: 

Ebin = Tra IH — . (21) 

V av' / 

The velocity dependence is fit by the function 

1^ ^ cJl + (^)^T^n [l + ; (22) 

Ebin L Vc2Vbin/ J L \ V J i 

the five parameters are listed in Table ^ (the fits should not 
be extrapolated to velocities much higher than shown in the 
figure). At low velocity Ecap/Ebin rises as ln(l/«) because 
the energy change C decreases exponentially with impact 
parameter. At high velocity Ecap/Ebin decreases as a power 
of V, which is clearer for the binaries with mi /m2 > 1. 
The velocity at the transition between the logarithmic and 
power-law behavior, approximately C2Vbin, depends on the 
mass ratio in the same way as w (eq. |l^) . 

For most velocities and mass ratios the cross sections 
in Figure ^ agree to within a factor of two with the approx- 
imate cross sections of Pineault and Duquet (1993). There 
are some larger differences at high velocity for an equal- 
mass binary, and at v/Vhin — 0.01-0.1 when m\/m2 3> 1. 
The (?;/Vbin)"* behavior seen when mi /m2 ^ 1 also agrees 
with that found by Pineault and Duquet. 

The capture cross section rises with the binary eccen- 
tricity, but the dependence is weak. The difference in Ecap 
for circular and highly-eccentric binaries is only 20-30%, 
too small to matter for most applications. 

Panels (c) and (d) of Figure ^ show the cross section 
for a captured orbit to survive for at least A'^ revolutions, 
i.e. for the radius to pass through at least A*' -I- 1 minima 
before the star is expelled. The results for binaries with 
mi/m2 > 1 are fit weU by Ecap(A)/Ecap ~ A"^/^. Ev- 
erhart (1976) noticed this N~^^'^ scaling in the survival 
of comets scattered by the Sun- Jupiter system, and inter- 
preted it as resulting from a random-walk in the comet's 
energy, as in the gambler's ruin problem from probability 
theory (see Yabushita 1979, Quinn, Tremaine, and Duncan 
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Figure 2. The cross-section Scap for a circular binary to capture an orbit with initial velocity v, plotted for five mass ratios on (a) 
linear and (b) logarithmic scales. The solid lines in (a) and (b) are the fits to eq. (^2[); the dashed line varies as (f/Vbin)~*- Panels 
(c) and (d) show the cross section for an orbit to be captured for at least A'^ revolutions for the same five mass ratios, plotted in (c) 
at the lowest v for each mass ratio, and in (d) a,t v = w. The dotted lines in (c) and (d) vary as N~^^^. 



1990 for further discussion). The N~^^'^ scahng does not 
work as well if the binary is not hard or if mi ~ m2. 

The cross sections in Figure ^ place no limit on the 
apocenter of the captured orbit. Some of the stars con- 
tributing to Scap are captured into weakly-bound orbits 
with apocenters many orders of magnitude larger than the 
binary's semimajor axis. In a real galaxy those orbits will be 
perturbed by passing stars and the galactic potential before 
they return to the binary. But that should not change the 
hardening rate much because the contribution from weakly- 
bound captures is small, even when mi ^ m2. 

3.1.3 Close-encounter cross section 

The cross section for close encounters with the binary mem- 
bers is needed for applications to real problems where the 
bodies are not point masses. For a massive BH binary we 
need it to compute the rate at which stars are tidally dis- 
rupted by the BHs, and to estimate how those disruptions 
might change the hardening rate. 

Figure ^(a) shows the cross section E for a star to 
approach within a distance < r of either of the BHs, for a 
hard, circular binary with m\/m2 = 64. The cross section 
is plotted for two sets of experiments: in the first the stars 
were allowed to encounter the binary only once, even if they 



were captured; the second allowed as many encounters as 
necessary for the stars to be expelled. 

The cross section for the larger BH scales as E/Ebin ~ 
r for the single-encounter experiments because of gravi- 
tational focusing. For the multiple-encounter experiments 
E is larger but the increase is important only for r/a < 
m-ilmx. The reason the increase is unimportant for r/a> 
m-ilm\ is that when a captured star orbits a binary with 
mi 3> m2 the star's distance of closest approach to mi re- 
mains nearly constant while its energy undergoes (approx- 
imately) a random walk. This is well known in cometary 
dynamics, where comets diffuse in energy at nearly con- 
stant perihelion (e.g. Duncan, Quinn, and Tremaine 1987). 

The cross section for a close encounter with the smaller 
BH is different. For the single-encounter experiments grav- 
itational focusing is important for r/a < ra2/m\ but not 
for r/a> m2/mi, so E/Ebin scales as r or as depending 
on whether r/a is smaller or larger than m^/mx. For the 
multiple-encounter experiments the cross section is larger 
for all values of r, not just for r/a < m2/mi. 

The distance r/a = m2/mi has a special importance 
for m2 if mi ^ m2 because the velocity of a star orbit- 
ing m2 at that distance equals Vbin. Figure |^(b) shows the 
cross section E for such encounters as a function of the 
mass ratio. For the single-encounter experiments E/Ebin ~ 
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Figure 3. Close-encounter cross sections for a hard, circular binary {v equals the lowest value in Fig. |l| for each mass ratio). Panel 
(a) shows the cross section for an orbit to approach within < r of m\ or m2 . The lines bounding the upper shaded region are for mi : 
the lower line results when the orbits encounter the binary only once; the upper when they have as many encounters as necessary for 
them to leave with positive energy. The lower shaded region has the same meaning, but for m2. The dotted lines vary as r"^ and r~^; 
the dashed line is at r/a = m2/m\. Panel (b) shows the cross section for an orbit to approach within r < {m2/m\)a of m2, for both 
single encounters (open circles) and multiple encounters (filled circles). Panels (c) and (d) show the differential hardening rates with 
respect to the distances of closest approach to mi and m2; the five lines are for mi/m2 = 1 (dotted), 4, 16, 64, and 256 (dashed). 



(m2/mi)^; for the multiple-encounter experiments cap- 
tures raise that scaling almost to E/Ebin ~ m^jmx. 

Panels (c) and (d) of Figure ^ show the differential 
hardening rates with respect to the distances of closest ap- 
proach to mi and m2, normalized so that the area under 
the curves is unity. The largest contribution to the hard- 
ening comes from orbits that pass both BHs at a distance 
not much smaller than the semimajor axis. When mi ^ m2 
there is a wide tail in the left of panel (d) , but the contribu- 
tion from close encounters with m2 is still a small fraction 
of the total hardening rate. 

In a real galaxy there will be two complications that 
can change these results. If weakly-bound captured stars 
are perturbed by nearby stars or the galactic potential they 
will not return to the binary in such a way as to keep 
their distance of closest approach to mi nearly constant. 
That would increase the difference between the single- and 
multiple-encounter cross sections for mi. But if the cap- 
tured stars are perturbed too much they might not return 



at all, which would reduce the cross sections for both mi 
and m2. The two complications tend to cancel for mi. 



3.H Distribution of final velocities 

The final velocity is the velocity of a star at infinity after it 
has been expelled by the binary. We need their distribution 
to compute the mass ejection rate. 

Everhart's (1968, 1969) work on the scattering of 
comets by planetary systems is relevant to the distributions 
to be considered here. Everhart used an approximate conic- 
matching procedure to derive the probability h{U) dU for 
the energy change U = AS* to lie in the interval dU after 
a single encounter between the comet and the planet. The 
distribution has three parts, which Everhart called A, B, 
and C. Parts A and B are for the small and intermediate 
energy changes and are fit well by (A) a Gaussian and (B) 
h{U) ~ 1/\U\^. Part C is for the large energy changes re- 
sulting from rare, close encounters with the planet. In parts 
A and B, h{U) depends only on \U\, but that symmetry is 
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broken in part C where energy gains are more frequent than 
energy losses. 

Panel (a) of Figure ^ shows the final- velocity distribu- 
tion for a hard binary from scattering experiments done in 
a manner similar to Everhart's, so that "final" means after 
a single encounter with the bina ry. If a star was captured 
its final velocity was set to — y'2|_Bt|, with Et measured 
when the star began returning to the binary for a second 
encounter; otherwise the final velocity was set to \/2El at 
the end of the integration. The figure shows the cross sec- 
tions E for the final velocity to be greater than or less 
than —Vf for some positive Vf . 

For the binaries with m\/m2 ^ 1 there is a range of 
velocities for which E is symmetric (depends only on |t;f|) 
and varies as E ~ !/«/■ This corresponds to Everhart's 
part B. The hardening rate would be nearly zero for these 
binaries if multiple encounters were not allowed because 
the positive and negative contributions would nearly can- 
cel. The symmetry is not as good for binaries with equal 
masses for which the star is more likely to gain energy. That 
is why the N~^^'^ scaling in Figure did not work so well 
for equal-mass binaries. The asymmetry that Everhart pre- 
dicted for part C is not clear in the figure, perhaps because 
the statistics are poor when mi ^ for the rare, close 
encounters with the mass m2 (the mi/m2 = 256 results 
come from lO'^ orbits, but that is still not enough). 

Panel (b) shows the final-velocity distribution when 
the stars are allowed to encounter the binary as many times 
as necessary for them to be expelled. The E ~ scal- 
ing from the single-encounter experiments is raised to ap- 
proximately E ~ ^/^fy but the probability of a star being 
expelled with — Vbin is still small if mi 3> m2. 

Panel (c) shows the differential hardening rate with 
respect to the final velocity, normalized so that the area 
under the curves is unity. The velocity at the maximum 
scales with the mass ratio in the same way as w (eq. |l^) 
and is approximately 1.75w. There is a wide tail to the 
right of the maximum if mi 3> m2 but the high velocities 
contribute little to the hardening. In fact the hardening 
rate for mi ^ m2 can be computed quite accurately by 
considering just the positive velocities from panel (a) and 
multiplying the result by two, i.e. by assuming that the 
captured orbits eventually get expelled with the same dis- 
tribution of final velocities as for the orbits that are not 
captured (this works only for very hard binaries). 



3.1.5 Discussion 

We can now explain why a hard binary hardens at a con- 
stant rate that is independent of its mass ratio. The ex- 
planation given by Roos (1981) is incorrect. It implies that 
the dominant contribution to the hardening for a binary 
with mi 3> m2 comes from orbits that have close encoun- 
ters with m2 and are expelled with high velocity. But the 
scattering experiments show that the dominant contribu- 
tion comes from orbits that do not have close encounters 
with the BHs and that are expelled with a velocity Vf ■r^ w. 

Consider a typical orbit that starts with a low velocity 
V, passes at a distance r ~ a from the two BHs, and leaves 



with a gain to its kinetic energy. The energy gain results 
mainly from the interaction with the smaller BH if mi 
m2 because the larger BH acts as a fixed potential. The 
interaction force of magnitude F ~ Gm2/a^ acts for a time 
At ~ (a^/GMi2)^^^ to produce a velocity change Av ~ 
FAt ~ (m2/Mi2)Vbin and a corresponding energy change 
A^;. ~ V^nAv ~ (m2/Mi2)Vfin- That gives C ~ m2/n ~ 
1, which is sufficient to give a hardening rate Hi with no 
dependence on the hardness and almost no dependence on 
the mass ratio. 

If this same derivation is repeated for a high-velocity 
star {v > Vbin) it gives a hardening rate that rises as 
Hi ~ instead of falling as Hi ~ as it should. 

That is because the derivation ignores the orbits that lose 
energy, which tend to cancel the ones that gain energy. 
The cancellation removes four powers of v and is the rea- 
son the hardening rate is so difficult to measure at large 
V by the Monte Carlo method. For a hard binary there is 
no cancellation because the orbits that lose energy in the 
first encounter are captured and eventually expelled with 
an energy gain. It is not surprising then that the hard/not- 
hard transition occurs at the velocity w where the binary 
begins capturing stars effectively. 



3.2 Mass ejection 

To measure the ejection rate we need an ejection criterion, 
i.e. a velocity Uoj such that a star with initial velocity v is 
counted as ejected if it is expelled with final velocity Vf > 
Ucj. The conventional escape-velocity choice, Vcj = 2%/3o", 
leads to a problem for a Maxwellian distribution because 
0.7% of the stars have initial velocities v > Uoj and will 
be counted as ejected if they receive any energy from the 
binary, no matter how little. We therefore choose «ej = 
max{1.5w, 2\/3cr}; the results do not depend sensitively on 
the numbers 1.5 and 2\/3. Let Fej{x,v, a) be the fraction 
of stars incident upon the binary with impact parameter x 
and initial velocity v that satisfy this criterion. The ejection 
rate is then 

J= — / dv Anv'^ f{v, a) — 4n / dx xFaj{x,v,a). (23) 
H Jo " Jo 

The integral over the velocity distribution is evaluated nu- 
merically after the inner integral is determined from final- 
velocity distributions like those shown in Figure ^. 

The ejection rate is plotted as a function of cr/ Vbin 
for five mass ratios in Figure ^ At low velocity J rises 
as ln(l/(T); at high velocity J falls, first gradually, then 
precipitously. The velocity dependence is fit by the function 



j2Vbir 



In 



(24) 



the five parameters are listed in Table q. The parameters 
give a close fit to the data in the figure but are erratic. Note 
that the velocity at the bend in the curves in panel (b) is 
not fit well by j2Vbin- 

The ejection rate for a hard binary can be estimated 
by noting that close encounters with the binary give a mean 
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Figure 4. Distribution of final velocities for a hard, circular binary (v equals the lowest value in Fig. ^ for each mass ratio). S is the 
cross section for an orbit to leave with velocity > Vf (solid lines), or to remain bound with energy < —v^jl (dotted lines). In (a) the 
orbits encounter the binary only once; in (b) as many times as necessary for them to leave with positive energy. Panel (c) shows the 
difi^erential hardening rate for the experiments in (b). The mass ratios •m\/m2 in (a) and (b) are as shown in (c), increasing from right 
to left. 
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0.5284 


0.8108 



Table 3. Parameters for fits to J (eq. 241) for a circular binary. 



energy change of (C) ~ 1. It then follows from the defini- 
tion (H) of C that a binary must interact with about its 
own mass in stars to shrink by a factor of e. But "interact 
with" does not mean the same as "eject." A binary that is 
not hard interacts with many stars but ejects few of them. 
And even a hard binary need not eject its own mass to 
shrink by a factor of e if it gives some stars much more 
energy than others. 

Figure (^) shows that a hard binary with mi/m2 3> 1 
expels few stars with Vf > Vbin- So why does J not decrease 
rapidly with mi/m2 in the left half of Figure (|^? Because 
although the fraction of stars expelled with Vi > Vbin does 
decrease rapidly with mi/m2, the fraction expelled with 
Vf > w does not, and it is that fraction that determines the 
ejection rate when the binary first becomes hard. 



3.3 Eccentricity growth 

Ki is more difficult to measure than H\ because of the can- 
cellation that occurs in the numerator of equation ( p^ ) . The 
B — C distribution is wide and nearly centered on the ori- 
gin with a mean {B — C) that is 10-100 times smaller than 
the deviation about the mean. Consequently we know much 
less about eccentricity growth than we know about harden- 
ing. Roos (1981) tried to measure K\ with only 500 orbits 
per measurement. He found K\ — 0.2 ± 0.2 for a hard bi- 
nary with e = 0.6 and concluded that the eccentricity could 
increase. Mikkola and Valtonen (1992) used 10'* orbits per 
measurement to study the dependence of Ki on the eccen- 
tricity and hardness of an equal-mass binary. Their results 
are accurate for hard binaries, for which they found pos- 
itive growth rates with a maximum of Ki = 0.19 ± 0.04, 
but have large error bars for v ~ Vbin- 

The results derived here use 10 - 100 times more or- 
bits per measurement than Mikkola and Valtonen used (10^ 
per measurement at low stellar velocity, 10*' at high veloc- 
ity) and use quasi-random numbers rather than random 
numbers to further reduce the statistical errors. The large 
number of orbits required has two practical consequences. 
First, only a small number of velocities and mass ratios can 
be examined. Second, the results should be applied only to 
problems where m* is much smaller than mi and m2, for 



Massive Black Hole Binaries 11 




10-3 ^ 



0.2 
0.15 
i< 0.1 
0.05 




_ 1 1 1 1 




1 1 1 1 1 1 _ 


: nCOl 


oQ.l 3^ 




r« 0.032 


• 0.32 




: * 1.0 
















0.2 0.4 0.6 0. 
Eccentricity e 



Figure 5. Mass ejection rate J for a circular binary, plotted on 
(a) linear and (b) logarithmic scales. The five lines are for mass 
ratios (increasing from right to left) mi/m2 = 1, 4, 16, 64, and 
256. 



Figure 6. Eccentricity growth rate for five initial stellar veloc- 
ities V, for binary mass ratios mi/7712 = 1 (a) and 16 (b). The 
lines are the fits to eq. (^). 
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Table 4. Parameters for fits to Ki (eq. |25[ ). 

otherwise the mean behavior can get lost in the dispersion 
about the mean (as often happens in N-body experiments). 

The measurements of K\ are plotted in Figure |^ for 
five initial velocities and two mass ratios. The eccentricity 
dependence for each choice of v and m\/m2 is fit by the 
function 

K^{e) = e{l-e^f° (ki + k2e); (25) 

the three parameters are listed in Table ^. K\ must ap- 
proach zero in the limits e = and e = 1: the former be- 
cause of the conserved Jacobi constant; the latter because 
of the (1 — e^) in equation (pl|). 

Consider first the results for an equal-mass binary. At 
large stellar velocity K\ is small and is negative for some 
eccentricities. As the velocity is lowered Ki rises and be- 
comes positive for all eccentricities. By u/Vbin = 0.032, Ki 



has converged to its limit for a hard binary. That limit gives 
a maximum growth rate of Ki ~ 0.2 near e = 0.7, consis- 
tent with the results of Mikkola and Valtonen (1992). The 
main difference between these results and theirs is at high 
velocity, where they did not identify negative values for K\ . 
Those values can reduce a binary's eccentricity when it first 
starts to harden. 

The results for m\/m2 — 16 differ from those for 
an equal-mass binary in two ways: the maximum Ki is 
about 40% larger, and the hard-binary limit is reached at a 
lower stellar velocity. In later applications we will need the 
growth rate for other mass ratios. We get that by assum- 
ing that the velocity w (eq. |l^ determines the transition 
to the hard-binary limit and that this is the only impor- 
tant dependence on the mass ratio, i.e. by interpolating the 
results for m\/m2 = 16 at the appropriate value of v/w. 

The derivation of the eccentricity growth rate for a 
Maxwellian distribution is more cumbersome than it was 
for the hardening rate; it is given in an appendix. The re- 
sulting K is the same as Ki at low stellar velocity but 
does not fall to zero as fast as K\ at high velocity. The 
appendix also derives K from Chandrasekhar's dynamical- 
friction formula and shows that it greatly overestimates the 
true growth rate. 



4 APPLICATION OF RESULTS TO MASSIVE 
BLACK HOLE BINARIES 
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4.1 Introduction 

We consider a galaxy core with uniform density p and ve- 
locity dispersion a. The core mass and radius, if needed, 
are computed from 



47r 3 



9a" 
47rGp 



1/2 



(26) 



Model I, for a large galaxy, has p = 10"^, a = 300, Mc = 
7.7 X 10®, and rc = 120 (dimensional quantities are given 
in units of Mq, pc, yr, and km/s). Model II, for a small, 
high-density galaxy, has p = lO'', cr = 100, Mc = 9.0 x 10*^, 
and To = 1.3. The BHs are assumed to enter the core and 
form a binary with initial eccentricity eo and semimajor 
axis ao, the latter chosen so that Vbin = cr/5. That is when 
the integration starts. 

The equations for da/dt and de/dt for three-body scat- 
tering are combined with those of Peters (1964) for gravi- 
tational radiation: 
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(27) 
(28) 



The equation for dMaj/dt is integrated to give the ejected 
mass but is ignored by the other equations, i.e. the density 
and velocity dispersion are held fixed. 



4.2 Characteristic length and time scales 

It is helpful to first consider how the length and time 
scales for the evolution vary with p, a, and the BH masses 
(see also BBR). The binary forms at a separation Ob = 
rc(Mi2/Mc)'^''^ where the enclosed stellar mass equals M12. 
It does not become hard until w > y/Sa, which happens at 
about 



6*7712 
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(29) 



From then the binary hardens in the time 
th = 



GpaH 



■ 4.3 X 10 yr | 



(}^\ (^) ri^ 



(30) 



where 77 ~ 16 is a typi cal hardening rate for a hard binary 
(recall that H is ■\/2/n times smaller than Hi). This should 
be compared with the time for a circular binary to merge 
through the emission of gravitational radiation: 
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(31) 



The two are equal at the semimajor axis 
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(32) 



The binary orbital velocity at a^r is close to the geometric 
mean of the velocity dispersion and the speed of light: 
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(33) 



This tells us the hardness at agr , what we previously called 
o-/Vbin: 



1/2 



(34) 



Two things should be noted about the length scales. 
First, flgr is large enough that tidal disruptions cannot 
change the hardening rate by much unless the BHs are 
very small. The disruption radius about the larger BH for 
a star like our Sun is 



which gives rt/flgr — 4.4 x 10~* for the numbers in equa- 



tion ( |32[ ) . That is too small to matter because the main con- 
tribution to hardening comes from orbits that do not have 
close encounters with the BH. Accretion of the disrupted 
stars will not matter either (although accretion from some 
larger gas reservoir can make a difference when the BHs 
are close). Tidal disruptions by the smaller BH can sup- 
press the high-velocity ejections if m\/m2 3> 1, but those 
ejections do not contribute much to the hardening. 

The second thing to note is the ratio Oh/agr, which 
determines how many e-foldings the binary has to harden 
between the time it becomes hard and the time radiation 
takes over: 



(103) (le) 



(36) 



The larger ah/ugr, the more mass the binary ejects and 
the more its eccentricity grows. For this example ah/flgr ~ 
exp(3.8) if mi — m2, but the ratio is smaller if mi is small 
or if m2 <^ mi . 



4.3 Evolution in a fixed galaxy 

The equations have been integrated for a number of mass 
combinations in the two galaxy models. Figure ^ shows the 
hardening time and the ejected mass for an initial eccen- 
tricity of eo — 0.1 (the eccentricity does not grow by much 
for this choice of eo). 

The hardening proceeds through three stages: the first 
ends when the binary becomes hard at a = a^, the second 
starts there and ends when gravitational radiation takes 
over at a = agr, and the third is the final merger stage. 
The gradual increase in the hardening time during the first 
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stage is because of the log term in equation ([L9|). All mass 
combinations harden along the same diagonal line during 
the second stage because all hard binaries harden at the 
same rate. The elapsed time equals the hardening time on 
that diagonal. The merger time in the final stage increases 
as mi decreases. The separation between the stages varies 
with the masses as expected: the smaller mi, or the smaller 
m2/mi, the less room the binary has to harden between ah 
and Sgr. Some binaries reach Ogr before they become hard. 

The hardening curves in panels (a) and (b) are sim- 
ilar to that shown in Figure 1 of BBR. The main dif- 
ferences are that BBR used a constant hardening time 
from Chandrasekhar's dynamical-friction formula during 
the first stage, and assumed that all binaries become hard 
at Vbin = cr. The hardening time used here must match onto 
Chandrasekhar's time at a large separation, but not until 
the log term in equation (|l^) matches the usual Coulomb 
logarithm. 

The ejected mass is never more than a few times M12 
and is less if the BHs are small or if mi/m2 3> 1. It can 
be estimated by Jgi- ln(ah/agr)/2 if ah/ttgr 3> 1. The mass 
ejected by a binary with mi/m2 ^ 1 can be much larger 
than m2 even though it is not for an equal-mass binary; so 
a BH of mass m ejects more mass if it merges with A'' BHs 
of mass m/N than if it merges with another BH of mass 
m. The mass displaced from the center of a real galaxy 
that starts with a density cusp will be larger than the Mcj 
computed here, because some is displaced when the BHs 
first approach the center and some is nearly but not quite 
ejected. 

The dependence on the initial eccentricity is shown in 
Figure ^ for one of the binaries from Figure 0. If eo < 0.3 
the eccentricity hardly grows before gravitational radiation 
drives it to zero; if eo ~ 0.1 the eccentricity goes down from 
the start. If eo > 0.3 the eccentricity grows, but not by 
much. This conclusion, which Mikkola and Valtonen (1992) 
also reached for an equal-mass binary, is true for all the 
binaries in Figure Q The merger time is reduced if eo is 
large, by about a factor of 20 if eo ~ 0.9 for the binary 
in Figure ^ The reduction could be larger if the harden- 
ing stalls during the late stages. Makino et al. (1993) say 
their N-body experiments show that massive BH binaries 
form with large eccentricities, but those experiments use 
unrealistic galaxy models and atypical initial conditions. 
Polnarev and Rees (1994) give arguments for why the ini- 
tial eccentricity should be small. 

4.4 Changes for a galaxy that is not fixed 

We can estimate when mass ejection is likely to affect the 
evolution by considering the mass of stars in the unper- 
turbed galaxy that pass within a distance r of a point mass 
M12 at the center: 

Mh<,.).IO«.(^)^(l + ^^). (3.) 

As the binary hardens, A/cj rises and M{rp < a) falls. If 
those two masses become equal the binary will have ejected 
nearly all the stars that can have a close encounter with it; 
further hardening must then wait for new stars to diffuse 
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Figure 8. Eccentricity evolution. The solid line is the same as 
in Fig. I^a); the other lines assume initial eccentricities of 0.3 
(dotted), 0.5, 0.7, and 0.9 (dashod-dotted). 

back into the "loss-cone" orbits by two-body relaxation. 
The point where Afoj = M{rp < a), marked by the open 
circles in Figure ^ typically occurs before the binary has 
ejected even one tenth of the mass that it has to eject to 
merge. The hardening rate will be reduced by mass ejec- 
tion before this point is reached. That will not change the 
Mcj(a) and e(a) relations but it will change the time scale 
for a{t), perhaps making the merger time longer than the 
age of the galaxy for some of the binaries in Figure ^. 

Mass ejection does not bring the hardening to a com- 
plete stop even in the absence of two-body relaxation be- 
cause the binary does not remain fixed at r = 0. A single 
particle of mass Afi2 would wander from the center of the 
unperturbed galaxy with an amplitude 

....;]^.«.01pe(^)(ig^)"". ,38, 

The mass scaling suggests that wandering will be more im- 
portant for small BHs, but a binary with large BHs can 
wander too once it ejects stars from the center because 
there is then no restoring force to keep it fixed. The impor- 
tance of wandering for the hardening rate is best studied by 
N-body experiments; all we can say here is that the merger 
times in Figure |^ are undoubtably too short. 

In summary, our study of the restricted three-body 
problem has answered some of the questions posed at the 
start. It has given a complete description of how the hard- 
ening, mass ejection, and eccentricity growth depend on the 
properties of a massive BH binary and a uniform galaxy 
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Figure 7. Hardening time \a/a\ and ejected mass Mej versus semimajor axis for BH binaries in two galaxy models; in model I (panels 
a— d), p = 10^ and c = 300; in model II (panels e and f), p = 10^ and a = 100. The initial eccentricity is 0.1. In (a) and (b), 
mi = m2 = 10* (solid line), 10^ (dotted), and 10^ (dashed). In (c) and (d), mi = 10* and m\/m2 = 1 (soUd), 4, 16, 64, and 256 
(dashed-dotted) . Panels (e) and (f ) are similar to (c) and (d) but use galaxy model II with mi = 10^ . The filled and open circles mark 
the points where a = ah (filled) and where the evolution could stall because of loss-cone depletion (open) . 



core in which it is embedded. It has cleared up some 
confusion resulting from mistaken applications of Chan- 
drasckhar's dynamical-friction formula. And it has allowed 
us to study binary evolution in uniform galaxy cores with 
some simplifying assumptions. But it has left two big ques- 
tions that affect the merger time: what the initial eccen- 
tricity is and by how much loss-cone depletion reduces the 
hardening rate. It has not allowed us to study binary evo- 
lution in realistic galaxy models with density cusps. And it 
has given only a crude estimate of the changes induced in 
such models by the evolution. These questions will be the 
subject of paper II. 
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APPENDIX A: ECCENTRICITY GROWTH 
RATE FOR A MAXWELLIAN DISTRIBUTION 

The derivation of K from Ki is more cumbersome than the 
derivation of H from Hi because Ki is a function of two 
variables (e and v), because we do not have as much infor- 
mation on its velocity dependence, and because the energy 
and angular momentum changes in equation ( [l^ must be 
averaged separately. We proceed as follows. We ignore the 
eccentricity dependence of Ix{C), which we know is small, 
and use formula ( |l6| ) for the velocity dependence. The value 
of Ix{B) is found at any eccentricity and velocity from the 
assumed Ix{C) and from linear interpolation of the mea- 
sured Ki{e,v). The quantities Ix{C)/v and Ix{B)/v are 
averaged over the Maxwellian distribution and substituted 
into equation ( p^ to give the growth rate K. There is some 
arbitrariness at high velocity, where we set Ki{e,v) = if 
V > 2Vbin, but that does not matter for our applications. 

The growth rate was computed in this way at 13 values 
of (r/Vbin and fit by the function (p^). The parameters for 
the fits are plotted in panels (a), (b), and (c) of Figure Al; 
the resulting fits are plotted in panels (d) and (e) for the 
two mass ratios. K a nd K i are about the same for a hard 
binary; there is no \plpi^ difference as there was for H. K 



does not fall to zero as fast as K\ at large stellar velocity, 
and negative values are less prominent for K. 

These results can be compared with the prediction 
from Chandrasekhar's dynamical-friction formula. The fric- 
tional force on a massive particle M moving with velocity 
V = ^aX is (eq. 7-17 of Binney and Tremaine 1987) 



dt 



erf(X) - ^ exp(-X2) 



(Al) 



Consider an equal-mass binary in units where G — Mi 2 = 
a — \. We ignore the factor of 47rG^AfplnA because it 
cancels in the ratio of the energy and angular momentum 
changes. Those changes are found by integrating over the 
unperturbed orbit (we use V here for the relative velocity 
of the two BHs): 



AS 



oc / dt2V'^F{V), 



dt F{V), (A2) 



where 



X 



erf (X/2) - — exp {-X^/4) 



(A3) 
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Figure Al. Eccentricity growtli rate K for a Maxwellian distribution. The top three panels show the fitting parameters (eq.p5|) 
for m\/m2 = 1 (filled squares) and 16 (open squares). The bottom three panels show the fit for mi/m2 = 1 (d) and 16 (e), and 
the prediction (for mi/m2 = 1) from Chandrasekhar's dynamical-friction formula (f); the six lines are for (from top to bottom) 
"■/Vbin = 0.002 5. .01. 0.04, 0.16, 0.64, and 2.56 (the top two lines are not distinguishable in (f)). The dotted line in (f) is the cubic 
approximation (K4). 



The integrals are evaluated numerically and substituted 
into equation ( p^ to give the growth rate -K'ch- Perturba- 
tion theory shows that for a hard binary 



le^ + 0{e')] (X»l), (A4) 



8 

which serves as a check on the numerical integration. 

The resulting growrth rate is plotted in panel (f). It 
differs from the correct growth rate in four ways: it falls to 
zero too fast as a rises; it rises to the hard-binary limit too 
fast as a falls; it approaches zero too slow as e approaches 
unity; and its maximum value for a hard binary is about five 
times too large. The last three differences together would 
greatly reduce the merger time for massive BH binaries if 
Kcii were the correct growth rate. But it isn't. 



